Discontinuous Molecular Dynamics for Rigid Bodies: Applications 
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Event-driven molecular dynamics simulations are carried out on two rigid body systems which 
diflFer in the symmetry of their molecular mass distributions. First, simulations of methane in 
which the molecules interact via discontinuous potentials are compared with simulations in which 
the molecules interact through standard continuous Lennard- Jones potentials. It is shown that 
under similar conditions of temperature and pressure, the rigid discontinuous molecular dynamics 
method reproduces the essential dynamical and structural features found in continuous-potential 
simulations at both gas and liquid densities. Moreover, the discontinuous molecular dynamics 
approach is demonstrated to be between 2 to 100 times more efficient than the standard molecular 
dynamics method depending on the specific conditions of the simulation. The rigid discontinuous 
molecular dynamics method is also applied to a discontinuous-potential model of a liquid composed 
of rigid benzene molecules, and equilibrium and dynamical properties are shown to be in qualitative 
agreement with more detailed continuous-potential models of benzene. Qualitative differences in 
the dynamics of the two models are related to the relatively crude treatment of variations in the 
repulsive interactions as one benzene molecule rotates by another. 



I. INTRODUCTION 

Although event-driven simulation studies of systems 
of molecules interacting via discontinuous potentials are 
common^, 0, IS S IE 0; most work to date has consid- 
ered fully-flexible models of molecular systems in which 
standard "bonded" potentials for the motion of bonds, 
bond angles, and internal dihedral angles are replaced by 
stepped, discontinuous versions of these interactions 
However, much of the internal motion of a molecule 
is largely unimportant in understanding many dynam- 
ical properties of condensed phase systems on long time 
scales. As the efiiciency of an event-driven simulation 
depends critically on the number of events to process 
per unit of simulation time, a great deal of computa- 
tional effort in simulations of large molecular systems 
is devoted to treating essentially unimportant internal 
motion of molecules, potentially limiting the scope of 
the simulation. Similar issues are relevant in more stan- 
dard, continuous-potential models of molecular systems. 
In this context, methods of incorporating partial inter- 
nal constraints in the dynamics of molecules for perform- 
ing efficient simulations of rigid systems with quaternions 
and matrix propagation schemes have been developed to 
extend the time scale or system size accessible to simu- 
lation. 

In the preceding paper P| , a general framework to carry 
out event-driven simulations of molecular systems inter- 
acting via discontinuous potentials in the presence of con- 
straints was outlined. This formalism is generally ap- 
plicable to both semi-flexible and fully-rigid bodies. In 
the case of a rigid molecular system, it was shown that 
the analytical solution of the free motion of arbitrary 



rigid bodies can be combined with efficient numerical al- 
gorithms for the search for event times and the use of 
appropriate collision rules to carry out event-driven or 
discontinuous molecular dynamics (DMD) simulations. 

In the present work, we apply this approach to two 
molecular systems of practical interest, methane and 
benzene. The lack of similarity in the mass distribu- 
tion of methane and benzene molecules, classified as 
spherical and symmetric top rotors, respectively, im- 
plies that the rotational motion in these two systems is 
quite different. First, a discontinuous-potential model 
of rigid molecular methane is constructed by analogy 
with a standard continuous-potential model based on 
Lennard- Jones interactions. Structural and dynamical 
characteristics of DMD simulations obtained at gas and 
liquid-like densities are then compared to their analogues 
in the continuous-potential system, and the computa- 
tional efficiency of the DMD approach is contrasted with 
standard molecular dynamics (MD) simulations. Second, 
a discontinuous model of a liquid of rigid molecules of 
benzene is designed to compare to a united-atom model 
of rigid molecules interacting with continuous potentials. 
Equilibrium and dynamical properties of both models are 
compared under similar conditions of density and tem- 
perature. The dynamics of orientational degrees of free- 
dom is examined and the few discrepancies observed in 
details of the angular motion are explained. 

The paper is organized as follows: The methane sys- 
tem is studied in Sec.^ starting with a discussion of gen- 
eral considerations pertinent for performing DMD as well 
as continuous-potential simulations of a symmetric-top, 
rigid molecular system in part A, followed by the presen- 
tation of gas phase results, the liquid phase simulation 
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results and analysis of computational efficiency in parts 
B, C and D, respectively. In Sec. 11111 simulations of a 
united-atom model of rigid benzene are presented. Con- 
tinuous and discontinuous potential models are presented 
in part A, followed by a discussion of the simulation de- 
tails, and the equilibrium properties results in parts B 
and C, respectively. A detailed analysis of the dynamical 
results, including a careful examination of the discrepan- 
cies associated with the differences between the models, 
is provided in part D of this section. Finally, conclusions 
are discussed in Sec. IIVI 



II. RIGID METHANE 
A. Specification of the model 

Consider a system of rigid methane molecules in which 
the four hydrogen atoms are held fixed with bond length 
d in a tetrahedral arrangement around the central car- 
bon atom. If the mass of the hydrogen atoms is taken 
to be the mass of one proton, m^, and the mass of the 
central carbon is 12 rup, the principal moments of iner- 
tia of the rigid rotor system are all equal and given by 
h = h ~ h = 8mp(P/3. As discussed in Sec. II of the 
preceding article^, statistical averages of a dynamical 
variable X{rM,PN), depending on the phase space coor- 
dinate {vmjPn), for the constrained system are defined 
by 

{X{rN,pN)) = ^ J dpNdrN\\Z{rN)\\ 
X X{rN,PN) p{rN,PN) 

xY\_(ii<7a{rN))S{&airN,PN)), (1) 

a 

where p{rj^,p^) is the probability density for the uncon- 
strained system, (Talr^) = are a set of 3n— 6 constraint 
conditions required to hold all distances between the n 
atoms in the molecule fixed, and Z is the partition func- 
tion for the constrained system, given by 

Z = J dpNdrN\\Z{rN)\\ p{rN,PN) 

X ]J(5(o-Q(rAr))5 {(^a{rN,PN)) ■ 

a 

Furthermore, in Eq. (Q, the matrix Z(r7v) is given by: 

1 daairiy) dapirN) 



dr.. 



(2) 



For this system, there are 9 total constraint conditions 
for each molecule and they can be chosen to be: 
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where the carbon atom is numbered 1, numbers 2 to 5 
refer to to hydrogen atoms, and a = 2(i-\/6/3 is the dis- 
tance between vertices in the tetrahedral structure. It is 
evident that the matrix Z of Eq. ||2J) is a 9 x 9 square 
matrix whose elements depend on the fixed distances 
between atoms in the molecule, as well as angles between 
them via dot products such as ri2 • ri^. Since the evolu- 
tion of the system holds these values fixed, the Z matrix 
is constant, and is invertible. Given that the Z matrix is 
constant, its determinant acts as a constant multiplica- 
tive factor for both the integral and the normalization Z 
in Eq. (Q), and hence factors out. 

For this rigid molecular model of methane, discontin- 
uous interaction potentials can be constructed by com- 
bining repulsive hard-core and attractive square-well in- 
teractions as follows. The repulsive hard-core interac- 
tion potential between site i on molecule a and site j on 
molecule b is given by 



ifr^<d,, 



if 



(3) 



where we have used the notation rfj — \rf—rj \ as the dis- 
tance between site i of molecule a and site j of molecule 6, 
and dij is the interaction distance. We will use the natu- 
ral choice d23 = ^24 = C?25 = ^34 = ^35 = ^45 = 0?HH 

(i.e., all hydrogen atoms interact the same way with 
one another) and di2 = dis — di4 = di5 = dcH (i-e., 
all hydrogen atoms interact the same way with carbon 
atoms), and we will also denote dn = dec- Choosing 
also dcH = ^(^cc + '^hh), the values of dec and c?hh can 
be interpreted as the diameter of the carbon and hydro- 
gen atom, respectively. However, a more general case is 
possible where the atom's effective diameter is dependent 
on the identity of its interaction partner. 

In addition to the hard-core repulsive interaction po- 
tentials above, we include an attractive square well (SW) 
between carbon atoms in the model of the form 



V-attr 
11 




(4) 



where the values of Vsw and dsw are adjustable parame- 
ters and the obvious condition dsw > dec must hold for 
this potential to have any effect. The attractive poten- 
tial allows for a condensation transition in a dense fluid 
of methane molecules at low temperature. 

The DMD results for this methane model presented 
below were obtained with the following parameter values: 
dec = 3.3 A, dHH = 2.75 A, dcH = 3.025 A, dgw = 5.1 A 
and Vsw = 1-824 kJ/mol. 

Given the form of the potentials in the model, the 
exact trajectory can be computed from arbitrary initial 
conditions. [1^ At any time t, the orientation of a given 
methane molecule is specified by a rotation matrix A(t), 
known as the attitude matrix, while its center of mass 
undergoes free linear motion in between impulsive events 
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that alter momenta in a discontinuous fashion. The atti- 
tude matrix and its transpose (t) can be used to trans- 
form between the lab and principal axis (body) frames 
and to calculate the Cartesian positions r,; (t) of an atom 
i in a molecule using 

= R{0) + V{0)t + A\t)-f„ (5) 

where fi is the constant position vector in the body 
frame, and R and V are the Cartesian vectors in the 
laboratory frame for the center of mass position and ve- 
locity, respectively. 

As discussed in detail in Rcf. 8, the time dependence 
of the attitude matrix is given by 

A(t)=P(i).A(0). (6) 



where P(i) is a rotation matrix itself which 'propagates' 
the orientation A(0) to the orientation at time t. For a 
spherical rotor in which all the principal moments of iner- 
tia are equal, the propagator matrix P{t) is particularly 
simple, 



P{t) - U(-(I;<), (7) 



where u) is the angular velocity vector in the body frame. 
In Eq. O, XJ{—ujt) — XJ{6h) is a rotation matrix rotat- 
ing an arbitrary vector by an angle 9 = —1^1^ around an 
axis n = <^/|<^| with components (ni, 712,713) according 
to 



"-1 + ('^2 + '^3) '^'^^ ^ ni7t2 (1 — cos 9) — 713 sin 9 7ii?i3 (1 — cos 9) + n2 sin t 
U(0n) = I nin2 (1 — cosS) -I- n3sin6' n2 + (nf + nfj cos 9 712713 (1 — cos 0) — ni sin | . (8) 
^773771 (1 — cos 0) — 7i2 sin 61 713712 (1 — cos 61) -I- Til sin 6* 7i| -I- (71^ -I- 712) cos 



For a spherical rotor, it follows from the Euler equations 
that a; is a constant vector determined by initial condi- 
tions, and hence the matrix U{—u3t) can be viewed as 
a rotation around a fixed axis by an angle that is a lin- 
ear function of time. Note that the moment of inertia 
tensor in the laboratory frame l{t) and in the principal 
axis frame I are related by l{t) = A^{t) ■!■ A{t), with 
I = diag(/i, /i, /i) for a spherical rotor, implying that 
l{t) is diagonal, constant, and equal to the moment of 
inertia tensor in the principal axis frame I. 

The angular momentum vector iD, defining the axis 
of rotation in P{t), remains constant until an impulsive 
force acts on it during a collision event. The times tc 
at which the impulses act arc determined by zeros of a 
collision indicator function fij{tc) = rfj* — dij, describing 
the time at which the distance rf^ between atom i of body 
a and atom j of body b attains a value at which there is 
a discontinuity in the potential energy (see Eq. Q). As 
explained in Ref. '3, the collision times can be computed 
using numerical grid-search methods using Eq. 

As discussed in Ref. H, the effect of the impulses on 
the Cartesian coordinates or angular velocities can be 
computed cither from a constrained variable approach, 
using the Z matrix, or using a rigid-body approach. In 
the rigid-body approach the impulse at collision time tc 
leads to a discontinuity in the center of mass momentum 
vector Pa and angular velocity vector of body a in 
the laboratory frame according to 



P' 

a 



Pa 



APa 
AuJn. 



(9) 
(10) 



where Pa = MVa is the center of mass momentum vector 
of body a, M is the total mass of body a, and P'a and 
(jj'a indicate the post-coUisional center of mass momentum 
and angular velocity vectors. Application of conservation 
laws leads to the following results: 



APa 



Au.a - -51-1 .(r-xr^^), 
where S is the magnitude of the impulse given by 



S = 



-6± V62-4aA$ 



(11) 
(12) 



(13) 



A$ is the height of the discontinuity in the potential at 
the interaction distance. 



1 



1 



AE''. -f AE''. 
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ab ^ ~ah 



and 



with 



AEt = n^.I-i.n^, 



In the equations above, all quantities are to be taken 
at time t = and r°i = rf — R'^ denotes the position 
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vector of atom i on molecule a relative to its center of 
mass il" and r°'j is the unit vector along the direction 
of the vector — rf connecting atom i on body a with 
its colliding partner j on body b. The physical solution 
of S corresponds to the positive (negative) root if 6 > 
(6 < 0), provided > 4a A$. If this latter condition is 
not met, there is not enough kinetic energy to overcome 
the discontinuous barrier and the system experiences a 
hard-core scattering, in which case A$ is replaced by 
zero and consequently S — —b/a. For the spherical rotor 
methane system, the fact that the moment of inertia ten- 
sor in the lab frame is proportional to the identity matrix 



leads to the simplification AE"; 



a.b 



With this model in hand, a DMD simulation can be 
carried out where four different types of events can take 
place, namely: hard-core collisions; square-well interac- 
tions; cell crossings; and virtual collisions associated with 
the truncation of numerical searches for collision events 
(described in detail in Ref. 0) . While the time at which 
cell crossings and carbon-carbon repulsive and attractive 
interaction events occur can be calculated analytically 
from the linearity of the center-of-mass motion, the times 
for the other atom-atom intcrmolecular interactions must 
be computed numerically using the grid-search method 
elaborated in Ref. ^. The use of cell subdivisions, local 
clocks, and a binary tree to manage the event calendar is 
a standard practice in this type of simulation and largely 
improves the simulation's efficiency 

For the purpose of comparison, a rigid molecular model 
of methane based on continuous interaction potentials 
was also implemented. In this model, all intermolecular 
interactions were assumed to be of Lennard-Jones form: 



12 
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where i?*; 



R 



■ ) 

mm/ 



and 



(14) 



Here, 



Ci and i?Jnin represent respectively the value and the po- 
sition of the minimum of the Lennard-Jones potential 
between atoms of the same kind, i.e. Viru). The param- 
eters Ci and are taken from Ref. where these 
and other parameters have been computed for atoms in 
a wide variety of molecules of biological interest and are 
intended to be used in standard condensed-phase simu- 
lations. The values of dij in the DMD model were in 
fact chosen to be comparable to i?J„i„. Standard values 
were used for all parameters in the continuous potential 
system: the C-H bond distance is 1.11 A, R^^^ = 3.60 A, 
= 3.13 A, i?™ = 2.66 A, ecc = 0.388 kJ/mol, 
ecH = 0.260 kJ/mol, and ehh = 0.176 kJ/mol. 

The continuous MD simulations were carried out in 
a standard fashion. The equations of motions were in- 
tegrated with a 4th order predictor-corrector algorithm 
with rotations represented by quaternion parameters. 
For efficiency, the forces and torques were calculated us- 
ing a cell structure where the index of molecules in each 



cell were stored in a linked list. In the simulation, the 
cell size was set to be equal to the cut-off value of the 
Lennard-Jones potential, taken to be 6.18 A. 



B. Low density results 

To study a low density methane system, DMD and 
continuous-potential simulations were performed for the 
models introduced above. In both cases, a system with 
512 molecules was simulated at a temperature T = 298 K 
(/cbT = 2.478 kJ/mol) and a density of 3.47-10-3 g/cm^. 




FIG. 1: From top to bottom, carbon-carbon, carbon-hydrogen 
and hydrogen-hydrogen RDFs for methane at a gaseous den- 
sity of 3.47 ■ 10"^ g/cm^ and a temperature of 298 K. The 
continuous lines represent the continuous MD model and the 
dashed lines represent the DMD results. The parameters of 
the model are given in the text. 
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AB I ri-r2 range I qab (DMD) | qab (MP) 

CC I 0-6.7 A I 0.18 ±0.01 I 0.18 ±0.01 

CH 0-6.7 A 0.71 ±0.01 0.74 ±0.01 

HH 0-6.7 A 0.71 ±0.01 0.73 ±0.01 



TABLE I: Number of nearby atoms associated with the var- 
ious peaks of the RDFs for the low-density methane fluid 
presented in Fig. The values are calculated using Eq. 1151 . 



Temperature is defined here as twice the kinetic energy 
divided by the number of degrees of freedom in the sys- 
tem. In both types of simulation, initial points in phase 
space were drawn from a so-called iso-kinetic ensemble 
where the kinetic energy is fixed according the temper- 
ature, while the potential energy can vary. In the DMD 
simulations, however, it turns out that the resulting fluc- 
tuations in the potential energy are so small (less than 
half of a percent) that a trajectory at one fixed energy 
may be used as representative of the whole iso-kinetic 
ensemble. 

In order to give correct time-correlation functions, the 
dynamics calculated in the simulations should in prin- 
ciple take place at constant total energy. However, the 
continuous-potential simulation exhibits a systematic en- 
ergy drift, which would violate the conservation of en- 
ergy. Therefore, it was supplied with a thermostat in the 
form of a periodic and infrequent (e.g. every picosecond) 
rescaling of the velocities. For consistency a similar in- 
frequent rescaling of the velocities could be performed 
in the DMD simulations, but this has a negligible effect 
because of the above mentioned small fluctuations of the 
potential energy in the DMD simulation. 

The DMD simulation, after equilibration, was run at 
constant energy of 7.37 kJ/mol. The pressure was cal- 
culated using the standard impulsive limit formula[T(Tj 
and was found to be 5.303 ± 0.009 bar. The continu- 
ous MD simulation was performed at constant temper- 
ature by using the above mentioned velocity-rescaling 
thermostat. The total energy and pressure values were 
7.236 ± 0.05 kJ/mol and 10.44 ± 0.09 bar, respectively. 

Although ideally, one would like to simulate a system 
at standard pressure (1 bar), the pressure in a simulation 
at constant volume and energy is very sensitive to small 
changes in the parameters of the interaction potential. 
This is true for both the DMD and the standard con- 
tinuous potential simulation technique. For this reason, 
pressures which do not differ by an order of magnitude 
(such as those above) may be considered as being in rea- 
sonable agreement. 

In Fig. n the intermolecular carbon-carbon, 
carbon-hydrogen, and hydrogen-hydrogen radial dis- 
tribution functions (RDFs) obtained from the DMD 
simulation are shown as a function of radial distance, 
and compared to those obtained from the continuous 
simulation. From the lack of correlation evident in these 
RDFS, it is clear that the structure of the fluid is that 



of a gas, which is the expected phase at this density 
and temperature. Note that the agreement between 
the carbon-hydrogen and hydrogen-hydrogen RDFs is 
quite good, while the carbon-carbon RDF exhibit a clear 
difference in the shape of the only peak that appears. 
The difference between the two models in the peak 
structure for the carbon-carbon RDF is to be expected 
since there is no interaction between molecules separated 
by a distance greater than the square-well parameter 
(5.1 A) in the DMD model. In the continuous model, on 
the other hand, the interaction distance extends up to 
6.18 A. 

A quantitative comparison of the RDFs can be carried 
out by calculating the average number qAB of atoms of 
type B that are close neighbors of an atom of type A, 
where qab is given by 

qAB = 47rpB / r^gAB{r)dr, (15) 

Jri 

where ps is the number density of type B and qab (^) is 
the radial distribution function for the pair AB. With 
the choice of ri — and r2 — 6.7 A, one finds that 
Ice — 0.18 from the carbon-carbon RDF for both the 
DMD and continuous-potential models, indicating that 
the average number of near neighbor molecules around a 
given molecule is essentially the same in the two models 
at low density. A similar quality of agreement of all qAB is 
observed for all other pair types, as is evident in Table 

C. High density results 

A more challenging test for the discontinuous model 
is in the high density and moderately low temperature 
regime where the system exhibits liquid behaviour and 
more complex structure. The conditions simulated here 
correspond to a density of 0.347 g/cm^ and a tem- 
perature of 126.8 K {k^T = 1.064 kJ/mol). Under 
these conditions, the average total energy in the DMD 
and continuous-potential model simulations was found 
to be — 5.583 kJ/mol and —6.11 ± 0.05kJ/mol, respec- 
tively, with corresponding pressures of 61 ± 10 bar and 
102 ± 58 bar (which are in reasonable agreement in the 
sense explained above). 

In Fig. 121 the intermolecular carbon-carbon, 
carbon-hydrogen and hydrogen-hydrogen RDFs are 
plotted versus radial distance. The correlation functions 
show a highly-ordered fluid and are consistent with 
a liquid phase. Note that most details in the RDFs 
obtained for the continuous model are well reproduced 
in the DMD model. 

For a quantitative comparison, we again examine the 
average number of near neighbors qAB, as defined in 
Eq. (|15|l . The results are presented in Table ITU From the 
values associated with the first peak in the carbon-carbon 
RDFs, we see that every molecule is surrounded on aver- 
age by about 13 other methane molecules in both models. 
Furthermore, the second peak integrates to about 46, a 
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value that is slightly larger than the value of 44 found 
with the continuous model. However, given that this dis- 
tance is well beyond the square-well interaction range, 
the agreement is quite good and suggests that excluded 
volume effects are the predominant factor in the struc- 
tural packing in this simple liquid system, while details 
of specific interactions are relatively unimportant. 

It also is worth noting that the agreement between the 
models found in the structural functions at both densities 
for the same temperature and similar pressure values sug- 
gests that they would have a similar phase diagram. An 
exhaustive analysis of the phase diagram of this model 




(A) 



FIG. 2; From top to bottom, carbon-carbon, 
carbon-hydrogen, and hydrogen-hydrogen RDFs for methane 
at a liquid-like density of 0.347 g/cm^ and a temperature of 
126.8 K. The continuous lines represent the continuous MD 
model and the dashed lines represent the DMD results. The 
model's parameters are given in the text. 



AB 


ri-r2 range 


q.AB (DMD) 


QAB (MD) 


CC 


0-6.19 A 


12.9±0.05 


13.4±0.05 


CC 


6.19-10.14 A 


45.7±0.05 


43.6±0.05 


CH 


0-6.33 A 


56.7±0.05 


53.7±0.05 


HH 


0-6.6 A 


58.0±0.05 


59.6±0.05 



TABLE II: Number of nearby atoms associated to the various 
peaks of the RDFs for the methane fluid presented in Fig. |5] 
The values are calculated using Eq. (1151 . 



is, however, beyond the scope of this paper. 

In addition to the structural properties of the DMD 
model, event-driven simulations also enable dynamical 
correlations to be computed in rigid-body systems. In 
fact, unlike continuous-potential systems, the DMD ap- 
proach allows dynamical properties to be calculated from 
essentially exact (i.e., up to machine precision) trajecto- 
ries due to the simplicity of the form of the interaction po- 
tential. Fig.Oshows the results obtained for the normal- 
ized time autocorrelation function of the center of mass 
velocity (VACF) of a given molecule in the fluid, C{t) = 
(V(0) • V{t))/{V{0) ■ V{0)), for the continuous and 
discontinuous potential models under liquid conditions. 
Considering the fundamental differences in the nature of 
motion in the two models, the agreement at a qualita- 
tive level of C(t) between the two models is somewhat 
surprising. Nonetheless, it is evident that the short-time 
behavior, in which the VACF in the DMD model decays 
more quickly than in the continuous-potential model, and 
the degree of anti-correlation in the velocities or "cage" 
effect at longer time scales differ appreciably. 

Because of the qualitative similarity in the VACFs 
for the two models, the value of the self-diffusion co- 
efficient, D = [kBT/m) C{t)dt, agrees reasonably 
weh: D is found to be (0.89 ± 0.03) • lO^^cm^/s and 
(3.10 ± 0.06) • lO^^cm^/s for the DMD and continuous 
MD model, respectively. The somewhat lower value for 




r(ps) 



FIG. 3: Velocity time correlation function for the case of liq- 
uid methane of Sec. Ill CI The continuous and dashed lines 
correspond to the MD and DMD results, respectively. 
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1000 10000 
Number of molecules 

FIG. 4: CPU time per picosecond of real time as a function of 
the system size for the low density methane case of Sec. urn 
The dots correspond to the DMD values whereas the open 
squares are the values obtained in a continuous MD simula- 
tion. The solid and dashed lines join the DMD and MD values 
and show the expected linear scaling with system size. 
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FIG. 5: CPU time per picosecond of real time as a func- 
tion of the system size for the high density methane case of 
Sec. Ill CI The dots correspond to the DMD results and the 
open squares correspond to the continuous MD simulation. 
The insert shows the ratio MD/DMD of the CPU time values 
as a function of system size. 



the diffusion constant in the DMD code might be due to 
the strong center-of-mass attraction in the DMD model, 
in contrast to the continuous model in which the attrac- 
tive potential is distributed over the atoms. The stronger 
attractive potential in the DMD model may cause the 
molecules to stay bonded for a longer time, thus decreas- 
ing the mobility. 



D. Efficiency 

Given the complexity of constructing an event-driven 
simulation for rigid systems, it is natural to wonder 
whether it is worthwhile to carry out such calculations. 
Such concerns can only be addressed by considering the 
relative efficiency of event-driven versus standard sim- 
ulations as a function of the system size and of the 
physical conditions of the simulation. In order to as- 
sess the relative computational efficiency of the DMD 
simulation, the CPU time needed to simulate one pi- 
cosecond of real time dynamics was assessed for the 
DMD and continuous-potential models for different sys- 
tem sizes and densities. Of course, such a comparison de- 
pends sensitively on a number of factors that have little 
to do with the methods themselves, such as the choice of 
programming language, computer architecture, compiler 
choice and degree of compiler optimization. In order to 
minimize such external effects, the DMD code was writ- 
ten in a high level language (C-I--I-) while the code for 
the continuous-potential simulations was adapted from 
Ref. ITfil and was written in C. Both codes were compiled 
with the same compiler, on the same computer cluster, 
and with the same non-specialized level of optimization. 
Since object-oriented code in C-|— I- is generally consid- 
ered to be less efficient than C because of the additional 
overhead of utilizing classes, it is likely that optimiza- 
tions in coding style and changing the coding language 



to C would enhance the relative efficiency of the DMD 
code reported here. 

It is also clear that the efficiency of a continuous MD 
simulation is largely dependent on the size of the time 
step chosen to integrate the equations of motion, i.e., 
the larger the time step the more rapidly the simula- 
tion propagates the system. However, a large integration 
time step leads to numerical instabilities in the MD tra- 
jectory, something that can only be tolerated with an 
extensive use of thermostats. The use of thermostated 
dynamics has its own limitations since it artificially bi- 
ases the "true" dynamics and influences the calculation 
of dynamical properties. The effect on the dynamics of 
thermostats is generally unknown, and hence artificial 
means of stabilizing trajectories must be handled with 
great care when looking at time dependent phenomena. 
In our simulations, the time step was chosen so that the 
effect of thermostatting the system was not evident in any 
of the dynamical correlation functions examined. This 
procedure resulted in a time step of 3 fs for the low den- 
sity runs and 0.9 fs for the high density simulations of the 
continuous-potential system. The efficiency of the simu- 
lations in the discontinuous model depends on the grid 
time interval used in the root search, and here a value 
of 20 fs was used for all simulations. In general, one can 
use the angular velocities of the colliding molecules to 
extract an average oscillation frequency of the collision 
indicator function and estimate the appropriate grid time 
interval for the simulation 12] from this information. In 
such an approach, the angular velocities, and hence the 
value of the grid interval, vary with temperature in much 
the same way that the optimum time step value changes 
in a continuous MD simulation with temperature. In- 
stead, we have chosen to use /cb?^ as the energy unit in 
the simulation so that the average angular velocities in 
the body frame have the same numerical value regardless 
of the temperature. This approach has the benefit that 
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FIG. 6: The relative speed-up due to the use of virtual colli- 
sions as measured by the ratio of the CPU time for simulating 
liquid methane for 1 ps without and with virtual collisions as 
a function of the number of time points in the search interval. 
The circles, triangles and squares correspond to the system 
with 1000, 3375 and 8000 methane molecules, respectively. 



the numerical value for the optimal grid time interval is 
also independent of the temperature. 

Figure 0] shows the CPU time in seconds required by 
the two methods to simulate 1 ps of real time dynamics 
for the low density system as a function of the number 
of molecules. Note that the the discontinuous method is 
faster than the MD simulation by two orders of magni- 
tude for all system sizes. This result is understandable 
on the basis that the dynamics in a low density system is 
dominated by free molecular motion which, as discussed 
earlier, can be integrated and propagated with large "fi- 
nite" time steps. This is to be contrasted with the or- 
dinary MD method in which an (ideally infinitesimally) 
small time step must be utilized in order to numerically 
integrate the dynamics and conserve energy. We also 
observe in Fig. ^ that both methods scale linearly with 
system size due to the use of cell lists and a cut off inter- 
action distance in the continuous model. 

Figure shows the CPU time needed by the DMD 
and MD methods to carry out 1 picosecond of real time 
dynamics for the denser system at 0.347 g/cm'^. The 
insert shows the CPU time ratio MD/DMD and indicates 
that, even at this density, the DMD method is more than 
2 times more efficient than the continuous MD method. 

One clear disadvantage of performing DMD simula- 
tions on rigid bodies as opposed to simple hard spheres 
is the necessity to find collision times using numeri- 
cal root-search methods. The evaluation of distances 
between possibly interacting atoms is computationally 
costly, and is often the most demanding task in a DMD 
simulation of rigid systems. To reduce the computational 
work of finding collision times, it is convenient to truncate 
the numerical search of collisions that extend beyond a 
certain maximum time interval, which will be called the 
search interval, since events far in the future are not likely 
to be executed. To implement a truncation scheme, col- 
lisions found in the search interval are scheduled as such. 



while a truncated search is scheduled as a virtual colli- 
sion (see Ref. 0). If the molecules for which the search 
for a collision event was truncated do not have an event 
over the previously searched interval, a virtual collision 
event is executed for the pair in which the root search 
is continued from where it was previously stopped. In 
this way, the numerical effort is primarily focused on the 
most likely events (i.e. the ones happening within the 
search interval) and the number of grid evaluations in 
which the Cartesian positions and the distances between 
all atoms in the pair of molecules are calculated, is di- 
minished. Even though the truncation of the numerical 
search for collision events reduces the computational load 
of this type of event, the introduction of new events in- 
creases the size and complexity of the binary tree used to 
manage events. As a result, the scheduling of new events, 
which typically scales as A^logA^, where N is the num- 
ber of events in the tree, becomes more demanding. Since 
the computational effort associated with the data man- 
agement increases while the effort of finding collisions 
decreases, the optimal size of the time interval for the 
virtual collisions is a function of the size of the tree, and 
hence the overall size of the system. If the search inter- 
val is very large, very few virtual collisions are scheduled 
and the efficiency of the simulation converges to that of a 
simulation that docs not include virtual collision events. 



Fig. shows the relative gain in efficiency due to the 
use of virtual collision events for the systems consisting 
of 1000, 3375 and 8000 molecules as a function of the 
maximum number of evaluations of the collision function 
in the search interval. Clearly, the use of virtual collision 
increases as this ratio is decreased. For the system with 
1000 molecules, an intensive use of virtual collisions (i.e. 
one or two grid evaluations) doubles the computational 
efficiency with respect to a simulation not utilizing these 
events. Smaller systems (not shown) exhibit the same 
trend and have similar relative efficiencies. For a very 
large system, however, the gain in efficiency is smaller 
due to the fact that a large binary tree slows down the 
processing of events. Indeed, for the system with 8000 
molecules, it is more efficient to limit the number of vir- 
tual collisions scheduled by extending the search up to 
two or three grid evaluations instead of performing a sin- 
gle grid evaluation per root search. The effect of the 
increased load in data management can be clearly seen 
as the system size increases. Nonetheless, the scheduling 
of truncated root searches leads to significant improve- 
ments in simulation efficiency for all system sizes exam- 
ined where the numerical search for collision times, as 
opposed to data management, represent the bulk of the 
computational work in the simulation. 
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III. BENZENE 
A. Molecular model 

The methane system considered in the previous sec- 
tion has particularly simple rotational dynamics due to 
the high degree of symmetry of the methane molecule. 
However the formalism estalalished in Ref. allows sim- 
ulations of rigid molecules with arbitrary mass distribu- 
tion to be performed. A simple example of slightly more 
complicated rotational dynamics is that observed in sym- 
metric top molecules, such as benzene, in which two of 
the principal moments are equal and the third is distinct. 

For a system of rigid benzene molecules, there are a 
total of 13 relevant sites in every molecule; the center 
of mass, six carbon sites, and six proton sites, where all 
atoms lie in a single plane with the carbon sites forming 
a closed ring. However using a united atom approach, 
in which every carbon-hydrogen pair is represented as a 
unique site, the number of relevant sites in every molecule 
is reduced to 7. If the radial distance from the center of 
mass to each united "carbon" atom in the plane of the 
molecule (taken to be the XY plane) is d and the mass 
of the united atoms is set to men = ISmp, then one 
finds that the principal moments of inertia are Ii = I2 = 
SmcH^^ = 39mpd'^ and I3 = 2/i. 

Typically, the intermolecular interactions between 
the united atoms in benzene are assumed to be of 
Lennard- Jones form, leading to an interaction en- 
ergy that depends strongly on the relative orientations 
of the interacting molecules. The simplest possible 
discontinuous-potential model possessing attractive but 
directional interactions consists of a combination of a 
spherically symmetric attractive well about each center 
of mass combined with purely repulsive interactions be- 
tween united carbon atoms. We therefore consider a 
united atom model of benzene in which the center of mass 
site, or site 1, is defined to have a repulsive hard- wall 
interaction with any other center of mass site, at the in- 
teraction distance dn while the united carbon sites repel 
one another at distance dec- The attractive interaction 
potential between the centers of mass of two molecules is 
defined as 

(-Vswi iirf^<dswi 

y.Ur ^^ab^ ^ ) ^y^^^ ^^^^ ^ ^ab < ^^^^ (^g) 

[ iirl'l>dsw2, 

where r^f is the distance between the center of mass 
of any two molecules a and b. This potential con- 
tains two attractive square- well regions and four parame- 
ters: two square- well distances dswi and dsw2, and two 
square- well depths ^514^1 and Vsw2- 

A DMD simulation with this model can be imple- 
mented in the same way as for the methane model in 
Sec. n, with a slightly more complex propagation matrix 
P(i) (see Eq. ©). For the symmetric top system, Pfi) 
can be written as the product of two rotation matrices[a|. 
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FIG. 7: From top to bottom, center of mass-center of mass 
(XX) and carbon-carbon RDFs for liquid benzene. The con- 
tinuous lines represent the results for continuous MD model 
and the dashed lines represent those for the DMD model (see 
Sec. inrHl for details). 



V{t)^V{-uo.ptz)-vl^^t^, (17) 

where ujp = [1 — /3//1) 0)2(0) is called the precession fre- 
quency, and Lq ^ I ■ uj[Q) is the initial angular momen- 
tum vector in the body frame. Unlike the case of the 
spherical rotor, the x and y components of the angular 
velocity in the body frame are time-dependent, and the 
moment of inertia tensor in the lab frame is no longer 
constant and diagonal. Thus, not only is the rotational 
motion more complicated, but the processing of collisions 
involves explicitly inverting the moment of inertia tensor 
to evaluate the change in angular velocities via Eq. H12|l . 

For the purposes of comparison, a continuous-potential 
united-atom model of benzene, in which the united-atoms 
interact via intermolecular Lennard- Jones interactions 
defined in Eq. 1^3}, was considered. Note that this model 
does not posses any explicit interactions between centers 
of mass and that the site-site interactions lead to an effec- 
tive orientational dependence in the interaction between 
molecules. 
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B. Simulation details 

Simulations were performed for both models in a box 
with 512 molecules at a density of 0.870 g/cm'^ with ini- 
tial conditions drawn from an iso-kinetic ensemble at a 
temperature of 298 K (fceT = 2.478 kJ/mol). The dis- 
tance d between carbon sites (or between the center of 
mass and every carbon site) was set to 1.48 A in both 
models of rigid benzene. 

The parameters in the DMD model used to obtain the 
results presented below are: dec — 3.25 A, dn — 2.0 A, 
dswi = 6.8 A, dsw2 = 12.8 A, Vswi = 1-982 kJ/mol 
and Vsw2 = 0.446 kJ/mol. The numerical search for col- 
lision events were carried out using a time grid of 19 fs 
and with only one time point in the search interval (i.e. 
extensive use of virtual collisions). The simulation was 
equilibrated with frequent rescaling of velocities to arrive 
at an initial condition drawn from an iso-kinetic ensem- 
ble, and later run with rescaling of the velocities only 
every picosecond. 

The Lennard-Jones parameters used for the contin- 
uous model simulations are = 0.457 kJ/mol and 
-^min ~ 3.695 A. The simulations were carried out in 
a standard manner using a predictor-corrector integrator 
on the translational and the orientational (quaternion) 
parameters. A time step of 1 fs and a cut off distance 
of 8.0 A for the truncation of the intermolecular inter- 
actions were used. With these parameters, the DMD 
method was 3 times more efficient computationally (mea- 
sured as the CPU time needed for the simulation of one 
picosecond of real-time dynamics) than the MD simula- 
tion. The continuous model was simulated at constant 
temperature with a thermostat to correct for the numer- 
ical energy drift. 

C. Static properties 

The average intermolecular energy values obtained for 
the DMD and the continuous model were —23.558 ± 
0.005 kJ/mol and -25.4 ±0.1 kJ/mol, respectively. The 
pressure values were 930 ± 25 bar and 360 ± 60 bar for 
the step potential model and the MD simulation, respec- 
tively. 

Fig. [7| shows the center of mass-center of mass and 
carbon-carbon radial distribution functions obtained for 
the DMD and continuous models of liquid benzene. The 
results agree with one another and with previous molecu- 
lar dynamics studiesfisj of liquid benzene at similar den- 
sities and temperatures. It is clear that the DMD model 
reproduces most of the structural details found with the 
continuous model in both radial distribution functions. 
In Fig.[7| the DMD result for the center of mass radial dis- 
tribution function shows two sharp decreases at distances 
of about 6.8 A and 13.8 A, which coincide naturally with 
the values of the square well attractive potential distances 
chosen for the model. Indeed, these distances were cho- 
sen such that the DMD model can match the continuous 
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FIG. 8: Normalized velocity time autocorrelation function 
for liquid benzene. The continuous line corresponds to the 
Lennard-Jones continuous potential and the dashed line cor- 
responds to the DMD result. 

radial distribution functions. 



D. Dynamical properties 

The normalized center of mass velocity time autocor- 
relation functions obtained for the DMD and MD models 
are shown in Fig. |S1 In spite of the differences between 
the models, the agreement between the correlation func- 
tions is very good. The values of the diffusion coefficients, 
obtained from integrating the functions in this figure, are 
(1.42 ±0.04) -lO-^cmVs and (1.62 ± 0.04) -lO-^cmVs 
for the DMD and MD simulations, respectively. The val- 
ues agree well, although both are smaller than the ex- 
perimental value of 2.27 •10~^cm^/s^2 reported at this 
temperature. 

The rotational motion in liquid benzene is more com- 
plex than the motion in liquid methane, and in other 
spherical rotor systems, due to the asymmetry of the 
benzene molecule. The asymmetry leads to non-trivial 
orientational components of the static structure in which 
the planes of nearby benzene molecules tend to be or- 
thogonal to one another. The asymmetry also leads to 
interesting correlations in rotational motion that can be 
examined by looking at normalized autocorrelation func- 
tions Cuj^. (t) of components of the angular velocity in the 
body frame, 

where the cuk are the k = x,y, z components of the an- 
gular velocity vector in body frame. In Fig. |51 the time 
autocorrelation functions of the x and z components of 
the angular velocity vector in the body frame are plotted 
versus time, where x and y are the principal axes in the 
plane of the molecule and z is the principal axis (Cg-axis) 
orthogonal to the plane of the molecule. Note that, due 
to symmetry, C^^{t) = C^^{t). 

As can be seen in Fig.|51 the x component of the angu- 
lar velocity time correlation function for the continuous 
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FIG. 9: Time correlation functions of the x and z components 
of the angular velocity vector in liquid benzene. The contin- 
uous line corresponds to the Lennard- Jones potential and the 
dashed line corresponds to the DMD result. The dotted lines 
are the results obtained with the hyperbolic tangent potential 
defined by Eq. llT^ . 



model exhibits a well defined minimum that is not repro- 
duced in the DMD model. Angular correlations of these 
components are short-lived in the dense system with a 
lifetime of less than a picosecond, due to the strong con- 
fining effect of molecular cages in the liquid. The min- 
imum in C^^ (t) calculated in the continuous potential 
model reflects an anti-correlation in the angular velocities 
due to molecules rebounding off of nearest neighbors and 
is reminiscent of anti-correlations observed in velocity au- 
tocorrelation functions. Interestingly, this phenomenon 
is not observed in the DMD simulation. 

In contrast, the autocorrelation function of the z com- 
ponent of the angular velocity has a much longer lifetime 
in both models, although C^^ (t) decays much too quickly 
in the DMD simulation. The longer lifetime arises due to 
the relatively free motion of rotations about the z-axis of 
a molecule, corresponding to a planar spinning motion. 

In order to understand the origin of some of the de- 
ficiencies of the discontinuous model, it is instructive 
to examine the interaction energy between two benzene 
molecules, lying in the same plane, as one molecule is 
rotated relative to the other around the z-axis for both 
models. This interaction energy for the continuous po- 
tential model is plotted in Fig. ^| In the discontinuous 
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FIG. 10: Repulsive potential energy (in units of ksT) for the 
continuous potential model between two benzene molecules as 
a function of their relative orientation measured by the rota- 
tion angle (in degrees) around the z-axis of one the molecules. 
The continuous, dashed, and dot-dashed lines correspond to 
coplanar benzene molecules at center-of-mass distances rn of 
6.0 A, 6.25 A, and 6.4 A, respectively. The circles correspond 
to molecules in orthogonal planes with rn = 6.0 A. 

model, the rotational motion is completely free unless 
the centers of mass of the molecules are close enough 
(rii < 6.25 A) so that a hard-core collision takes place. 
When the hard-core collision occurs, the z-component of 
the angular velocity in the body frame changes direc- 
tion and the rotation effectively reverses direction. In 
contrast, the repulsive energy in the continuous model 
increases up to a maximum where the distance between 
united-carbon atoms on the two molecules is minimized. 
The repulsive energy drops as the molecule rotates past 
this point, so that there is a "hindered" rotational in- 
teraction for the rotation of one molecule past another. 
Such oscillatory interactions, clearly absent in the DMD 
model, are responsible for the qualitative differences ob- 
served in Fig. 1^1 Note that the attractive term of the 
Lennard- Jones interaction can be effectively accounted 
for via center of mass interactions, and that its inclusion 
in this analysis will only further enhance the oscillatory 
structure of the interaction. An interaction analogous 
to the one shown in Fig. ^| can be included in the DMD 
model by adding finite repulsive shoulder interactions be- 
tween united carbon atom sites. 

Alternatively, the role of the hindered rotational in- 
teractions on the angular correlations can be isolated 
by systematic modifications of the continuous interac- 
tion potential to approach the impulsive potential limit 
and, therefore, the DMD model. For example, instead 
of using a Lennard- Jones interaction, repulsions between 
united carbon atoms can be modeled with a potential of 
the form, 

i^^^j) = "1 t^^li - rl')] . (19) 

The parameters ai and a2 can be independently assigned 
to control the height and the width of the repulsive poten- 
tial, respectively. The value of ai can be set large enough 
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FIG. 11: Time evolution of the z component of the velocity 
vector (top panel) and the x component of the angular ve- 
locity vector (bottom panel) of a benzene molecule involved 
in a collision with another benzene molecule. The solid lines 
correspond to the DMD simulation. The dot, dashed and 
dot-dashed lines correspond to the model interactions given 
by Eq. 1191 with a2 values of 5, 10 and 20, respectively. 



that the probabihty of any molecule going through the 
barrier is virtually zero. The hard-wall limit in which 
the interaction range approaches zero and interactions 
become instantaneous is formally obtained in the limit 
of infinite a2- Of course, such a potential cannot be inte- 
grated numerically, and comparisons of dynamics using 
numerical integration schemes and dynamics in the im- 
pulsive limit must use a large, but finite value of a2- 

As an example, consider the trajectories of two ben- 
zene molecules before and after their mutual interaction 
through the potential defined by Eq. H19I) . In Fig. llll the 
dynamical evolution of separate components of the linear 
and angular velocity vectors of one of the molecules in- 
volved in the collision is shown as a function of time for 
the DMD model. The trajectories in the impulsive limit 
are also included in this figure for comparison. Since the 
center of mass velocity vector, its components, and the 
z component of the angular velocity vector are constants 
of the motion, their behaviour is analogous to the one 
shown in the top panel of Fig. ^2 It is clear that the 
interaction time goes to zero and the final Vz value ap- 
proaches the DMD result as the value of 012 increases. 
The X and y components of the angular velocity vector. 



FIG. 12: Orientational time correlation functions Pi{t) and 
P2{t) in liquid benzene. The continuous line corresponds to 
the Lennard- Jones potential and the dashed line corresponds 
to the DMD result. 



on the other hand, follow an oscillatory dynamics and 
their time evolution after the interaction can be quite 
different than that observed in the impulsive limit, as 
can be seen in the bottom panel of this figure. Once 
again, one observes that the dynamics generated by the 
continuous potential approaches the DMD dynamics as 
the value of 02 increases. 

The time-correlation results obtained with continuous 
MD simulations of a system composed of 512 benzene 
molecules interacting via Eq. H19|) . at the same density 
and temperature as above, are also shown in Fig. |51 
These results are obtained with the values ai = 40 and 
a2 = 20 for the carbon-carbon repulsion. Not surpris- 
ingly, the new time correlation functions agree much bet- 
ter with results obtained with the DMD model than those 
of the Lennard- Jones model. This agreement confirms 
that the discrepancies observed in the angular velocity 
correlation functions in the discontinuous and continuous 
models are a consequence of the overly simplistic inter- 
action potential between united carbon sites. 

In addition to correlation functions of the angular ve- 
locities, one can also examine how long individual ben- 
zene molecules retain their orientation in a dense liq- 
uid. One common measure of orientational correlation 
involves computing ensemble averages of low order Leg- 
endre polynomials of the cosine of the angle between an 
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orientational vector of a molecule, such as the Cg axis 
(or z-axis), at two times. In Fig. ^1 the Pi{t) and P2{t) 
orientational time correlation functions, defined by 

Pi(i) = (i(0) •£(<)) (20) 

and 

P2it) = (^l{mo)-z{t)r'i}) (21) 

for the DMD and Lennard- Jones models are shown. 
The overall shape of these functions is in agreement 
with previous molecular dynamics simulations of liquid 
benzene|l5|. Although there is qualitative agreement in 
orientational correlations in the continuous and discon- 
tinuous models, again it is evident that orientational cor- 
relations typically decay more quickly in the discontinu- 
ous model than in the continuous model. It is likely the 
longer-lived orientational correlations in the continuous 
case arise due to the hindered rotational interactions that 
are poorly modeled in the simple discontinuous system. 

IV. CONCLUSIONS 

In this paper DMD simulations of two rigid molecu- 
lar systems, which differ in the symmetry of their mass 
distribution, were presented. The simulation methodol- 
ogy utilized is based on exact solutions of the free mo- 
tion of rigid bodies in combination with collision rules 
derived from conservation laws. In Sec. ^ a discontin- 
uous model of rigid methane was constructed in which 
methane molecules attract one another via a center of 
mass square- well potential, while the constituent atoms 
of different molecules repel each other through hard-core 
repulsions. DMD simulations of the spherical rotor sys- 
tem were carried out in low and high density regimes, and 
excellent agreement with standard molecular dynamics 
simulations was observed for all structural quantities. It 
was furthermore found that dynamical quantities, such 
as the self-diffusion coefficient, also agree well with those 
of standard MD simulations. 

In addition, the efficiency of DMD simulations relative 
to standard continuous-potential simulation methods of 
the rigid methane system was examined in order to as- 
sess the role of certain parameters such as the grid search 
size and the truncation interval for the root search in the 
DMD simulations. It was found that the DMD simu- 
lations were between 2.5 and 100 times more efficient 
than standard continuous MD simulations at densities of 
0.347 g/cm^ and 3.47 • 10"'^ g/cm^, respectively. Fur- 
thermore, it was found that it is generally most efficient 
to schedule virtual-collision events (i.e. truncated root 
searches) as frequently as possible for small systems and 
slightly less frequently for larger systems. The use of 
virtual-collision events typically leads to improvements 



in simulation efficiency of 25% for large systems and up 
to 250% for small systems. 

In Sec. mil a united-atom model with discontinuous 
potentials for rigid benzene was presented in which ben- 
zene molecules attract one another via center of mass 
interactions, while united carbon atoms repel one an- 
other via hard-core interactions. DMD simulations of 
the rigid benzene system at a liquid density were carried 
out using the simulation methodology appropriate for a 
symmetric-top rigid body. To the best of our knowledge, 
this is the first event driven simulation of a non-spherical 
rigid body, other than a linear dumbbell, where the dy- 
namics and collision rules are properly taken into account 
without approximation of the dynamics or collision rules. 
Static and dynamical quantities from the DMD simula- 
tion were compared with those obtained from a continu- 
ous potential benzene model, in which the united carbon 
atoms interact via a standard Lennard-Jones potential 
and found to be in good agreement. The quantitative 
discrepancies in dynamical correlations of angular veloc- 
ities and orientational directors between the models are 
readily explained in terms of the crude level of description 
in the discontinuous model of the hindered rotational mo- 
tion of one benzene molecule near another. Nonetheless, 
one can argue that the relevant physics is already present 
in the crude description of the interactions in the DMD 
model, such as long-lived orientational correlations and 
correlation times of lUz which are significantly larger than 
those of ujx and ujy . If one is interested in more detailed 
features of the dynamics, it is always possible to add more 
discontinuities in the interactions, at the expense of com- 
putational efficiency. In a sense, this is contrary to the 
spirit of the DMD method, which is aimed at efficiently 
simulating the main physical characteristics of a system. 

The major purpose of this work is to establish that 
event-driven simulations of rigid molecular systems can 
be carried out with impressive gains in efficiency over 
standard molecular dynamics methods. Of course, the 
magnitude of the gain in relative efficiency depends on 
a number of factors, such as the level of detail in con- 
structing a discontinuous potential model of the system, 
the system size, and physical parameters such as the den- 
sity. Nonetheless, our results indicate that the gain in ef- 
ficiency through the use of crude discontinuous potential 
models of condensed phase systems may allow larger sys- 
tems or longer time scales to be investigated than those 
currently accessible with standard simulation methods. 
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